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Abstract 

In this paper, a randomized numerical approach is used to obtain approximate solutions for a class of nonlinear Fredholm 
integral equations of the second kind. The proposed approach contains two steps: at first, we define a discretized form of 
the integral equation by quadrature formula methods and solution of this discretized form converges to the exact solution 
of the integral equation by considering some conditions on the kernel of the integral equation. And then we convert the 
problem to an optimal control problem by introducing an artificial control function. Following that, in the next step, 
solution of the discretized form is approximated by a kind of Monte Carlo (MC) random search algorithm. Finally, some 
examples are given to show the efficiency of the proposed approach. 
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Introduction 

It is well known that linear and nonlinear integral equations 
arise in many scientific fields such as the population dynamics, 
spread of epidemics, thermal control engineering, and semi- 
conductor devices. Integral equations are usually difficult to solve 
analytically. Therefore, numeric solutions of integral equations 
have been a subject of great interest of many researchers. In most 
cases, the problem of finding numeric solutions for the nonlinear 
Fredholm integral equations is more difficult than the problem of 
the linear Fredholm integral equations. The problem of solving 
numeric solutions for the linear Fredholm integral equations of the 
second kinds is one of the oldest problems in the applied 
mathematics literature. Many numerical computational methods 
including deterministic techniques and Monte Carlo methods are 
introduced in this area ([1-10]). These proposed methods of 
literatures [1-10] have been used to discuss the linear Fredholm 
integral equations. However, it is difficult to solve the nonlinear 
kinds with them, 



k(xj)F(<p(t))dt, a<x<b, 



(1) 



especially analytically, where g(x) is a known continuous function 
over the interval [a,b], F(-) is a known nonlinear function with 
respect to <P(x) and k{x,t) is a kernel function which is known and 
continuous too, at the same time bounded \k(s,t)\<M on the 
square D = {(x,t)\xe[a,b],te[a,b]} which M is the upper bound on 
the square D. Consequendy, our aim here is to find the unknown 
function <P(x) which is the solution of problem (1). 

Some deterministic techniques have been applied to obtaining 
solutions of the nonlinear integral equations during the past several 



years. Many different numerical and approximate techniques were 
introduced to obtain the solutions of nonlinear integral equations. 
Literature [1 1] applied a different method (DM) to solve the 
equation (1) by convert the problem to an optimal control problem 
based on introducing an artificial control function. The results of 
algorithm [1 1] depend explicitly on the selection of the starting 
interval which decreases the speed of convergence substantially. 
Discrete Adomian decomposition method (DADM) [12] arose 
when the quadrature rules are used to approximate the definite 
integrals which can not be computed analytically. Homotopy 
perturbation method (HPM) introduced by [13] was used to solve 
the nonlinear integral equations (1). Positive definite function 
method (PDFM) [14] is based on interpolation by radial basis 
functions (RBFs) to approximate the solutions of the nonlinear 
Fredholm integral equations. The accuracy of PDFM relies on the 
choices of radial basis functions. Triangular functions method 
(TFM) [15] was utilized as a basis in collocation method to reduce 
the solutions of nonlinear Fredholm integral equations to the 
solutions of algebraic equations by using the optimal coefficients. 
Harmonic wavelet method [16] was employed as basis functions in 
the collocation method towards approximate solutions of the 
Fredholm type integral equations. Optimal homotopy asymptotic 
method (OHAM) was introduced in literature [17] as a reliable 
and efficient technique for finding the solutions of integral 
equations. Literature [18] presented a computational method for 
solving nonlinear Fredholm integral equations of the second kind 
which is based on the use of Haar wavelets. Iterative method [19] 
was used to approximate the solutions of the nonlinear Fredholm 
integral equations (1) by transforming the integral equation into a 
discretized form. However, each of these methods has its inherent 
advantages and disadvantages. For some methods, such as iterative 
method, a good initial value must be given in advance. Otherwise, 
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Table 1. Algorithm 1. 




Step 


Operation 


1 


+ 1 T 

Choose t > 0, the maximum number TV of random search and step size b of random search. Select an initial solution y^ER" , where yg = (vocj'oi > ■ ■ ■ vVorc) 




and y 0 can be randomly selected by some probability distributions. 


2 


Set m = 0. 


3 


Compute *I>(yo). 


4 


If <t(yo)<(, go to step 10. 


5 


m = m + 1, if 777 — iV>0 then set b = b/2 and m = 0. 


6 


Generate 77+ 1-dimensional random variable yi = 0'io,}'n, ■ ■ ■ ,>'i H ) T from a normal distribution, that is y 1 ~N(y Q ,Z) l E = diag(/3,6, . . . 


7 


Compute fl^yi), if flKyitaflKyo). go to step 5. 


8 Set >„ = >! and O(y 0 ) = O(y 1 ). 


9 


Go to step 4. 


10 


If the stopping criterion is met, then stop 
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the obtaining results by the iterative method may diverge or 
converge to not true solutions of the equations. For some integral 
equations, the nonlinear function F(-) must satisfy convergence 
conditions of algorithm, otherwise interval of convergence would 
be local subinterval which is not global interval [19]. 

In this paper, we intend to present a random numerical scheme 
for obtaining approximate solutions of the nonlinear Fredholm 
integral equations (1) by random search algorithm and suppose 
that the discussed integral equations have one solution at least. At 
first, we transform the equation into a discretized form by a 
numerical method of integration, e.g. Simpson's rule. 

The present work is motivated by the desire to obtain numerical 
solutions of the linear Fredholm integral equations of second kind 
by using Monte Carlo method, which is introduced by Spanier 
and Gelbard [6], Veach and Guibas [7], Farnoosh and Ebrahimi 
[8], Zhimin et al. [9] and Doucet et al. [10]. But there have been 
only a few works that solved the nonlinear Fredholm integral 
equations such as (1) numerically by Monte Carlo methods. 

As a kind of numerical method, Monte Carlo method has a 
great many merits. For example, it has only a few moments that 
Monte Carlo method is limited to geometric restrictions, the 
convergence speed is irrelevant to the dimensions of the problem, 
the error is easy to determine, the program structure is simple, 
more flexible and easy to accomplish. Especially, when the solved 
nonlinear system has a higher dimension, deterministic technique 
is difficult to determine the solution of system. For numerical 



Table 2. Errors of numerical results <?„ 


of Example 1 and Example 2. 








Example 1: e„ 




Example 2: e„ 




n MC 


Symbolic solution 


MC 


Symbolic solution 


3 0.000317436 


0.000312506 


0.000779989 


0.000780082 


5 2.37532E-05 


2.33971 E-05 


8.86961 E-05 


8.85708E-05 


7 3.591 33E-06 


3.53349E-06 


2.081 52E-05 


2.081 37E-05 


9 8.69071 E-07 


7.67474E-07 


7.061 89E-06 


7.06065E-06 


13 1.14246E-07 


1 .07254E-07 


1 .46446E-06 


1.46353E-06 
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problems in a large number of dimensions, Monte Carlo methods 
are often more efficient than conventional numerical methods. 

The paper is divided into five sections. In section 2 we transform 
the equations (1) into a discretized form by Simpson quadrature. 
In section 3 Monte Carlo random search algorithm is introduced 
to solve the numeric solution of nonlinear algebraic system. In 
section 4 some numerical examples are solved by the proposed 
method. Section 5 concludes this paper with a brief summary. 

Simpson Quadrature Method 

We now approach the subject of numerical integration. The 
goal is to approximate the part of definite integral in the integral 
equations (1) over the interval [a,b] by evaluating integrand 
k(x,t)F{(p{t)) at some sample points. 

Consider a uniform partition of the closed interval [a,b] given by 
tj = a+jh of step length h = (b — d)/n,n>0 integer, and 
to = a,t„ = b, using Simpson quadrature, the integral term 
\ b a k(x,t)F((p(t))dt may be expanded as [20], 

n 

J2^k(x,tj)F(<P(tj)) + 0(h 4 ), (2) 

where Wj are weights and determined by the quadrature methods. 

Further more, the solution of integral equation (1) with 
quadrature methods is approximated by solving the following 
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Table 3. Errors <?„ 


of simulated results for methods PDFM [14], DM [IV 


and iterative method [19] 


in Example 1 


and Example 2. 




n 


Example 1: e„ 




Example 2: e„ 






MC 


[14]RBF(f,(x)=- n L T ) 


MC 


DM [11] 


[19](/i= 1/100) 


5 


2.4E-05 


1.7E-03 CAT =5) 


8.7E-05 






10 


8.7E-07 


2.1E-05 (AT =10) 


7.1 E-06 


2.0E-02 


7.0E-03 


15 


7.2E-10 


1.5E-08 (JV=15) 


7.6E-07 
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system, 



hold, then 



<P{xd = g( Xi ) + X ^ WjKxiJjWKtj)) + 0{h\ (3) 

7=0 

where xj = tjj = 0,1, • ■ • n. 

When n is odd, Simpson's 3/8 rule is applied to the last three 
subintervals and then Simpson's 1/3 rule is applied to the 
remained subintervals. 

Denoting the notation 

<p(Xi) =<P t ,g(xi) = gi, k(xi,tj) = kjj. 



If the truncation error of equation (3) is neglected, the nonlinear 
algebraic system can be obtained, as follows, 



P i = gi + X^2wjkijF0> J ). 
7=0 



(4) 



If let <Pi = <P{Xj) be an exact solution of equation (1) at point 
X = Xu then <P = {9q, ( P\, ■ ■ ■ ,<P n ) T is also a solution of equation (3). 
Denoting <f>* =(<P* 0 ,<P\, ■ ■ ■ ,<P*„) T as a solution vector of nonlinear 
system (4). 

The conditions that <P* approaches to exact solution <p will be 
given by the following corollary [19]. 

Corollary 1 Suppose the following conditions 



(i) — ^ v " exists for every te[a,b\, 
in) < c, \k(x,t)\ <M, MC< 



1 



W(b- a y 



lim ll^-^H, 

/i->0 



= 0. 



(5) 



where h = (b — a) /«, n > 0. 



Proof 

Please refer to literature [19]. 

For nonlinear algebraic system (4), Monte Carlo method will be 
used to obtain the approximate numeric solutions of (4) in the next 
section. The program structure of this method is simple, flexible 
and easy to accomplish. For problems whose required precision is 
relatively high, Monte Carlo method is a very good choice to 
achieve a good initial value. 

Monte Carlo Random Search Algorithm 

With regard to the nonlinear algebraic system (4), we introduce 
notation v, = <Pi,i = 0,1, . . . ,n that simplifies our work and system 
(4) can be rewritten as 



'/o6>o,J>i, 
/i6>o,j>i, 



■,y„)=o 
■,y„)=o 



(6) 



where fi(y 0 ,yi , • • ■ o>») = JV ~gi ~ * Y^j=o ' v k ' '' iy l 



Denoting y = (y 0 ,yi, 
objective function: 



' >JV) > by equation (6), we establish 



Table 4. Errors of numerical results e„ of Example 3 and Example 4. 





n 


Example 3: e„ 




Example 4: e„ 




MC 


Symbolic solution 


MC 


Symbolic solution 


3 


0.031940214 


0.031942343 


0.009557852 


0.009557697 


5 


0.004775551 


0.004748677 


0.00391024 


0.003910249 


7 


0.001225184 


0.00122517 


0.000631981 


0.000631981 


9 


0.000432635 


0.000432178 


9.44288E-05 


9.44288E-05 


13 


9.16911E-05 


9.16792E-05 


4.48755E-06 


4.48744E-06 
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Table 5. Errors <?„ of simulated results in Example 3 based on harmonic wavelet [16] and Haar wavelets [18] methods. 



n MC 


harmonic wavelets: (N,M) 


Haar wavelets: k =32 


8 1E-04 


9E-02 (N = 2,M = &) 




1 6 2E-05 


8E-03 (jV = 3,M=16) 


2E-02 


20 9E-06 


8E-07 (jV = 6,M=128) 
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Table 6. Errors e„ of simulated results for methods HPM [13], TFM [15] and OHAM [17] 


in Example 4. 




n MC 


HPM [13]((»(.v)a<I>5) TFM [15](m = 16) 


OHAM [17] 


4 8.3E-03 


16 9.2E-16 


1 .2E-06 2.0E-08 


2.5E-16 


20 2.8E-16 


doi:1 0.1 371 /journal.pone.01 03068.t006 



3>(y)=]T/; 2 (y). 



(7) 



Obviously, if y = (j'o,J ; i, • ■ ■ ,y„) is the solution of equation (6), 
then y also makes the objective function fl>(y) equal zero. Hence, 
the unknown solution y of nonlinear equation (6) can be 
determined in such a way that the objective function cD(y) is 
minimized. Namely, suppose that there is a vector 



,y*) such that 
<%*) 



min ®(y)<£ (e>0) 



(8) 



holds, then y* =(y^,y\, ■ ■ ■ ,y*) is not only an approximation 
solution of nonlinear equation (6), but also an approximation 
solution of nonlinear integral equation (1) at point 
X = Xj,(i = 0,l, . . . ,n) when the kernel function k(x,t) and nonlin- 
ear function F(<p(t)) satisfy the conditions of Corollary 1 . 

In this paper, numerical methods based on random sampling, 
especially Monte Carlo methods, can be applied to this 
optimization problem. 

For solving problem (6), we consider Monte Carlo random 
search algorithm in Table 1. 

Numerical Examples 

In this section, we are going to demonstrate some numerical 
results for solving <P(x) in the nonlinear problem (1). Therefore the 
following examples are considered and the numerical solutions are 



x 10 



-4 




Figure 1. Errors of MC random search with n=13 and variant value x in Example 5. 

doi:1 0.1 371 /journal.pone.01 03068.g001 
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Table 7. Errors <?„ of simulated results for Legendre wavelets method [21] in Example 5. 





// 


MC 


Legendre wavelets method 


Legendre wavelets method 


4 


7.6E-04 






16 


4.6E-06 


1.8E-03(M = 3,/c = 2) 


1.3E-04(Af = 4,/fc = 2) 


20 


1 .9E-06 
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Table 8. Error results for Example 6. 





n MC 


Method of [23] 


5 5.8E-03 


5.2E-04 


9 1 .9E-03 


1.3E-04 


17 5.5E-04 


3.2E-05 


doi:1 0.1 371 /journal.pone.01 03068.W08 



obtained based on Monte Carlo random search algorithm. In 
order to measure the efficiency of the proposed MC random 
search algorithm, the symbolic solution is also used to solve the 
following examples. 

Remark: In Table 2 and Table 3, results of error e„ of the 
symbolic solution here is referring to the truncation error of 
equation (3). 

Let <P(x) be exact solution of the nonlinear Fredholm integral 
equation (1), y* = {y\,y\, ■ ■ ■ ,y*) T (where y* =<P*,i = 0,\, ...,«) is 
an approximation solution obtained by using the given Table 1 
with a known £>0 and partition step size h and 

J'o = C)'oO'>'oi' ' ' ' >yon) T i s an m Mal value by random sampling. 
In all examples (except Example 3), we choose the initial value 
C/(0,1). In Example 3, we set up the initial value Vq~ (7(1,2). 
By reference [14], we define the maximum error to measure the 
efficiency of MC method for different values n, which is given as 



e„= max \<P{x t )- 

0<i<n 



(9) 



Example 1 

Consider the following nonlinear Fredholm integral equation 
[14] 



<P(x)-- 



1 



dt+x- 



, 0<x<l 



which has the exact solution <P(x) = x. 

Errors of the numerical results for different values n are given in 
Table 2. Table 3 lists comparison of results with PDFM [14]. 

In this example, for reference [14], as radial basis functions 
(RBF) take different forms, the accuracy based on RBF is different. 
For instance, as radial basis functions <^ 3 (x) = 1 /( 1 + x 2 ), the 
maximum error is 1.7 x ltr 3 with N=5. But for our proposed 
method in this paper, from Table 2, the maximum error is 
2.4 x 10~ 5 with the same value of n. 

Example 2 

Consider the following nonlinear integral equation [11,19] 



<P(x)-- 



(x + t)e m dt + ex - 1 , 0 < x < 1 . 




0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

X 

Figure 2. Approximation process of MC random search with n = S in Example 1 for solving <P(x). 
doi:1 0.1 371 /journal.pone.01 03068.g002 
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1.2 
I 

0. 

0.6£ 
f 

0.4 
0 

-0.2 



8 
8 



8 



Exact 
x Initial value 
o Approximation 



o.i 



0.2 



0.3 



0.4 



0.5 

X 



0.6 



0.7 



0.8 



Figure 3. Approximation process of MC random search with n = 5 in Example 2 for solving <p(x). 
doi:1 0.1 371 /journal.pone.01 03068.g003 



0.9 



For which the exact solution is <p{x) = x. 

One may find in Table 2 the comparison of exact and obtained 
approximate solutions for different values n. Comparison of results 
with DM [11] and iterative method [19] is displayed in Table 3. 

Comparison of the results of this example with those obtained in 
[1 1] and [19] shows the efficiency of the proposed MC algorithm 
more accurate. 

Example 3 [16, 18] 



yj) = (0,l/n,2/ft, • • • ,1) , symbolic solution is failure. But for 
Monte Carlo random search algorithm, the order of maximum 
error e„ is still 10 ~ 2 . 

Table 4 shows errors of the numerical results for Example 3 
with different values n. Table 5 shows that the proposed Monte 
Carlo method solves equation (1) more accurately than methods 
harmonic wavelets [16] and Haar wavelets [18]. 

Example 4 [13, 15, 17] 



<P(x)-- 



xt[<p(tfdt + e* - y ,J , 0 < x < 1 . 



<P(x)-- 



cos(7ix)sin(7i?)KO] sin(nx), 0<x<l. 



For which the exact solution is <P(x) = e x . In this example, for 



symbolic solution, if we set up the initial value = sin (nx) + - 



In this example, the analytical solution of this equation is 
20-v/39T 



-COS (rex) on interval [0,1]. 
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1.2 




' 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

X 

Figure 5. Approximation process of MC random search with « = 5 in Example 4 for solving <p(x). 
doi:1 0.1 371 /journal.pone.01 03068.g005 



Table 4 shows errors of the numerical results for Example 4 
with different values n. Table 6 describes comparison of maximum 
errors e„ with HPM [13], TFM [15] and ORAM [17]. From 
comparison of the results of Example 4, we see that MC method 
provides a better numerical solution for the nonlinear Freholm 
integral equation of second kind and is more preciser than 
deterministic methods HPM and TFM. The proposed MC 
method and OHAM [17] can have the same order of error, is 
10~ 16 as n is large. 

Example 5 [19] 



cp(x) = e x 



x+2 x + 2 



e x '<P 2 (t)dt, 0<X<1 



The error function e„(x) in Figure 1 shows the precision of the 
approximate solution in example 5. The result shows that the 
interval of convergence of the proposed MC method is global 
interval [0,1] and the order of error is 10~ 5 . But for reference [19], 
the interval of convergence of the deterministic method is local 
subinterval [0,0.01] and the order of error is 10 . 

Example 6 

Consider the following nonlinear integral equation [21] 



where g(x)=gi(x)-( 



cp(x) = g(x) + 
9 9 



xt(f> (t)dt, 



7 



1 



128 32e 16e 2 16e 4 



which has the exact solution <p(x) = e x . 
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u.o 

0.6 
0.4 
0. 



o 
8 



Exact 
x Initial value 
o Approximation 



x 

8 



-0.2 
-0.4 
-0.6 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

x 

Figure 7. Approximation process of MC random search with n = 4 in Example 6 for solving <p(x). 
doi:1 0.1 371 /journal. pone.0103068.g007 



8 



gi(x> 



*e[0,-), 



1 1 



The exact solution is <p(x) = g"i(x). Table 7 illustrates the 
numerical results of Example 6 using methods of [21] and the 
proposed MC. 



<p"(x)-e' p(x) =0, xe(0,l), <p(0) = <p(l) = 0. 
This problem has the unique solution 

<P(x) = - In (2) +2 1b (c- sec — ±- ), 



where c is a unique solution of C- sec — = and this nonlinear 



Example 7 

As the final example we consider a physical problem which is of two-point boundary value problem may be reformulated as the 



great interest in magnetohydrodynamics [22]. 

The proposed MC method is used to solve an integral equation 
reformulation of the nonlinear two-point boundary value problem 



integral equation 




* Exact 
x Initial value 
o Approximation 



8 



4 




0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

x 



Figure 8. Approximation process of MC random search with n=5 in Example 7 for solving <P(x). 
doi:10.1371/journal.pone.0103068.g008 
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<p( x ) = g ( x )+ f k(x,t)e^dt, xe[0,l], 
Jo 

where the free term g(x) = 0, kernel function 

k(x,t)=< x),x_t ^ ^ Green's function for the 

^ —x(\ — t),x<t 

homogeneous problem 

<p"(x) = 0, xe(0,l),H0)=«'(l) = 0. 



The errors results using the method of [23] together with the 
results obtained for errors e n by the proposed MC method are 
tabulated in Table 8. 

Approximation process of Monte Carlo random search 
algorithm of Examples 1-7 were drawn with given value of n in 
Figures 2-8. 

As was seen in Figures 2-8, numerical solution of equation (1) 
based on the proposed Monte Carlo random search algorithm is 
convergent in probability to the analysis solution and this proposed 
MC method is stable. 



Conclusion 

In this paper, a numerical method based on quadrature 
methods and Monte Carlo random search algorithm has been 
proposed to approximate the solutions of the nonlinear Fredholm 
integral equations. With this method, the problem of solving 
integral equation is reduced to a problem of solving a nonlinear 
system of algebraic equation. Illustrative examples are given to 
demonstrate the validity and accuracy of the proposed MC 
method. In Tables 1, 3, solutions of Examples 1-4 are presented 
for different values of n to show that the maximum errors e„ based 
on the proposed MC method and the direct solving method of 
equation (4) have the same order. The results obtained by MC 
random search algorithm are compared with solutions obtained by 
DM, DADM, HPM, PDFM, TFM, OHAM, harmonic wavelets, 
Haar wavelets, Legendre wavelets, iterative method and colloca- 
tion-type method. Tables 2, 4, 5, 6, 7 and Figure 1 show that the 
proposed MC method is reliable and efficient for finding the 
solutions of integral equations in Examples 1-7. Figures 2-8 shows 
that the propose MC method is convergent and stable. For 
concrete problems that the region of root of integral equation (1) is 
unknown, the proposed Monte Carlo random search algorithm is 
a good and reliable choice to obtain an initial value and achieve a 
high precision when we deal with the problems. By the way the 
program structure of this method is simple and easy to accomplish. 
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